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1  Introduction 

An  application  of  hybrid  optimization  techniques  to  airfoil  design  with  Navier-Stokes  equations  is  here  pre¬ 
sented  and  discussed.  The  greatest  obstacle  in  the  use  of  Navier-Stokes  equations  with  evolutionary  optimiza¬ 
tion  procedures  is  the  huge  computational  effort  required.  To  overcome  this  limit,  several  approaches  have  been 
suggested,  that  range  form  parallel  genetic  algorithms  [1,  2]  to  approximated  fitness  evaluation  [3,  4]  to  hybrid 
techniques  [5]. 

These  three  approaches  can  be  profitably  combined  and  an  example  will  be  here  given  related  to  transonic 
airfoil  design. 

The  lecture  is  organized  as  follows:  an  outline  of  hybrid  optimization  techniques  will  be  given  then  the  par¬ 
allel  hybrid  genetic  algorithm  will  be  presented.  Thereafter  an  application  example  will  be  reported,  and  finally 
some  ongoing  developments  about  the  hybrid  approach  based  on  gradient  computation  through  the  adjoint  to 
Navier-Stokes  equations  will  be  described. 


2  Hybrid  genetic  algorithms 

Hybridization  has  been  one  of  the  first  techniques  adopted  for  improving  genetic  algorithm  performance,  while 
keeping  the  desirable  flexibility  features  of  genetic  algorithms  [6], 

In  the  optimization  context,  hybridization  can  be  defined  as  the  mixing  of  two  or  more  search  techniques 
with,  possibly,  complementary  features. 

The  genetic  algorithm  itself,  even  in  its  simplest  implementation,  can  be  considered  as  a  hybrid  optimization 
technique,  where  selection,  mutation  and  crossover  cooperate  in  the  same  optimization  process. 

Figures  1  and  2  show  how  the  genetic  algorithm  can  be  extended  through  the  addition  of  a  hill  climbing 
operator.  In  the  first  scheme  the  hill  climber  receives  in  input  the  whole  population  belonging  to  the  current 
generation,  while  in  the  second  one  an  intermediate  selection  process  chooses  a  subset  of  elements  that  will  be 
refined  through  hill-climbing. 

A  hybrid  algorithm  requires  care  in  balancing  the  various  components  of  the  search  procedure.  Indeed,  if 
high-rated  solutions  arc  injected  in  the  population  at  an  early  evolution  stage,  this  may  adversely  affect  popu¬ 
lation  diversity  and  force  the  evolution  in  wrong  directions.  The  way  to  avoid  this  depends,  in  general,  on  the 
particular  optimization  technique  introduced.  If  a  gradient  based  technique  is  used  as  a  solution  refinement  op¬ 
erator,  it  may  be  useful  to  balance  the  improvement  rate  of  both  techniques  by  stopping  the  hill  climber  after 
a  few  iterations.  However,  the  iteration  number  choice  is  mostly  a  question  of  practice  and  experience,  and  it 
should  take  into  account  both  optimization  algorithm  and  problem  features. 

Another  point  that  requires  special  care  when  using  hybrid  techniques  is  encoding.  Often,  indeed,  a  special¬ 
ized  optimization  technique  relies  on  a  particular  encoding  of  the  problem  variables  that  may  play  a  key  role  in 
the  efficiency  of  the  method.  Therefore,  whenever  possible,  it  is  worthwhile  to  extend  to  the  genetic  algorithm 
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Figure  1 :  A  simple  hybridization  scheme  between  a  genetic  algorithm  and  a  hill  climber. 


Figure  2:  A  more  complex  genetic-hill  climber  coupling  scheme. 


the  encoding  of  the  specialized  technique.  On  the  other  hand,  the  simple  binary  encoding  used  by  many  genetic 
algorithms  may  cause  problems  to  the  specialized  algorithms  if  the  conversion  of  data  is  not  carefully  handled. 
For  example,  if  a  binary  genetic  algorithm  does  not  use  enough  bits  for  variable  encoding,  the  improvement 
obtained  using  a  gradient  based  hill  climber  may  be  lost  when  re-encoding  the  variables  in  binary  form. 

The  hill  climber  operator  adopted  here  is  a  gradient  search  based  on  BFGS  algorithm  [7].  Indeed,  the  fitness 
functions  of  the  presented  applications  arc  differentiable,  for  which  gradient  based  techniques  are  much  more 
efficient  to  locally  improve  a  given  solution.  The  genetic  algorithm  developed  adopts  a  bit  string  codification 
of  the  design  variables;  anyway,  this  does  not  prevent  the  use  of  operators  requiring  real  number  list  encoding, 
such  as  extended  intermediate  crossover  and  word  level  mutation  [8].  In  this  cases,  the  binary  string  is  decoded 
into  a  real  numbers  list,  the  operator  is  applied  and  the  modified  variables  set  is  encoded  back  into  a  bit  string. 
This  scheme  allows  the  use  of  a  free  mix  of  different  type  of  operators.  The  hybrid  genetic  algorithm  operates  as 
follows:  through  the  application  of  the  selection,  crossover  and  mutation  operators,  an  intermediate  generation 
is  created  from  the  current  one;  afterwards,  if  the  hybrid  option  is  activated,  some  randomly  chosen  individuals 
are  fed  into  the  BFGS  based  operator  to  be  improved,  and  then  introduced  into  the  new  generation. 
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3  Parallel  genetic  algorithm 

The  hybrid  genetic  algorithm  here  described  has  two  points  that  can  be  executed  in  parallel:  the  evaluation 
of  the  new  population  members  after  the  mutation  and  recombination  phase  and  the  evaluation  of  the  gradient 
through  finite  differences  in  the  gradient  based  hybrid  operator.  The  parallel  programming  model  adopted  relies 
on  shared  memory  multiprocessing  and  the  parallelism  is  implemented  at  the  thread  level  using  the  standard 
POSIX  thread  interface. 

The  same  code  base  works  on  a  SGI  POWER  CHALLENGE  system  with  16  R-10000  processors  and  on  a 
Linux  PC-cluster  with  eight  processors  and  the  MOSIX  clustering  software  [9]. 

MASTER  process 


Figure  3:  Parallel  evaluation  loop. 


The  algorithm  is  organized  following  the  master-slave  paradigm.  Figure  3  shows  the  architecture  of  the 
parallel  evaluation  loop. 

In  the  initialization  phase  that  precedes  the  first  execution  of  the  evaluation  loop,  the  master  process  creates 
a  pool  containing  a  number  of  processes,  equal  to  the  maximum  number  of  threads  available  for  the  computation 
(NPROC).  The  child  threads  arc  immediately  put  in  a  wait  state,  and  they  will  remain  in  such  a  state  until  they 
receive  a  “go  ahead”  signal  from  the  master  to  start  the  computation.  This  choice  avoids  the  inefficiency  of 
creating  a  child  process  every  time  a  computation  is  needed  and  killing  it  at  the  end. 
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When  the  master  process  enters  the  evaluation  loop,  it  splits  the  population,  of  size  NPOP,  in  slices  of  maxi¬ 
mum  NPOP/NPROC+ 1  elements  and  then  each  slice  is  assigned  to  one  of  the  child  threads.  Afterwards,  a  signal 
is  sent  to  the  child,  through  a  standard  POSIX  semaphore,  so  that  it  begin  to  evaluate  its  slice.  When  the  child 
terminates  its  computation,  it  sends  a  completion  signal  to  the  master  (using  another  POSIX  semaphore).  The 
master  waits  for  the  completion  of  all  child  processes  in  a  synchronization  point.  The  child  processes  arc  ter¬ 
minated  at  the  very  end  of  the  program,  when  all  the  evaluation  loops  related  to  each  generation  and  gradient 
evaluation  loop  have  been  completed. 

This  architecture  has  maximum  efficiency  when  each  thread  has  an  even  computational  charge.  If  this  is 
not  the  case,  the  computation  process  can  loose  efficiency  because  the  master  has  to  wait  at  the  synchronization 
point  for  the  completion  of  the  slowest  process.  A  partial  solution  to  this  problem  consists  in  activating  more 
threads  than  available  processors  (e.g.  8  threads  on  a  four  processor  machine).  If  this  is  not  possible,  then  the 
algorithm  should  be  modified  in  order  to  redistribute  dynamically  the  computational  charge  to  the  child  threads. 

4  Aerodynamic  flow  solvers 

Two  flow  solvers  with  different  level  of  accuracy  have  been  adopted  for  the  objective  function  computation  of 
the  presented  design  exercises.  The  first  one  is  the  in-home  developed  ZEN  Navier-Stokes  solver  [10]  with 
Baldwin  Lomax  turbulence  model  [11].  The  second  one  is  Drela’s  MSES  code  [12]  that  is  based  on  a  finite- 
volume  discretization  of  the  Euler  equations  on  a  streamlined  grid.  The  viscous  region  is  computed  using  an 
integral  boundary  layer  based  on  a  multi-layer  velocity  profile  representation.  The  inviscid  and  viscous  regions 
arc  coupled  using  displacement  thickness. 

The  boundary  layer  code  is  used  both  as  a  low-fidelity  solver,  and  as  a  helper  for  ZEN  code.  Indeed,  when 
constant  lift  coefficient  (q)  is  required,  a  suitable  angle  of  attack  for  ZEN  is  guessed  using  MSES;  after  a  first 
computation  with  ZEN,  the  angle  of  attack  is  corrected  using  the  dci/da  computed  with  MSES. 

5  RAE  2822  airfoil  optimization 

The  optimization  problem  requires  the  minimization  of  the  drag  coefficient  obj  =  cy,  with  control  on  lift  co¬ 
efficient  and  maximum  thickness.  The  design  point  is  fixed  at:  M  =  0.68,  Re  =  5.7  x  106,  q  =  0.56.  The 
starting  geometry  is  the  RAE  2822  airfoil.  Transition  is  fixed  to  x/c  =  0.01. 

The  airfoil  shape  is  defined  as  a  linear-  combination  of  an  initial  geometry  ijq  (x)  and  some  modification  func¬ 
tions,  fi(x),  i  =  1 ,n: 

n 

y(x)  =  y0(x)  +  ^2wifi(x)  (1) 

i= 1 


where  iq  are  the  design  variables. 

The  modification  functions  used  here  are  reported  in  Table  1  with  (  =  i/c  £  [0,1]. 

5.1  Optimization  using  simple  GA  and  Euler+BL 

A  first  optimization  run  has  been  performed  using  the  genetic  algorithm  without  hybrid  operators,  and  the  MSES 
solver  only.  The  geometry  was  represented  using  20  design  variables  (10  for  the  upper  surface  and  10  for  the 
lower)  chosen  among  the  polynomial,  Hicks-Henne  and  Wagner  functions.  Three  subsequent  runs  of  the  parallel 
GA  were  performed,  and  10  threads  were  active  in  each  run.  The  GA  parameters  are  reported  below: 


Variables  encoding 

32  bit 

Selection 

3  step  R.  W. 

Crossover 

one -point,  pc  =  1 

Mutation 

bit  level,  pm  =  0.01 

Population  size 

40 

Generations 

40 
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Hicks-Henne 

Legendre 

0.888  (1  -  0  v^e-13’28^ 

0.42  (l-O3^ 

0.57  (1  -  0 

0.946  (1  -£)3£ 

0.1  sin3  (tt£0-431) 

0.1  sin3  (7t£0'757) 

0.1  sin3  (V1-357) 

0.1  sin3  (tt£3-106) 

0.136  (1-03  (l2£  —  10£2) 

(1  -  £)3  (225£5  -  630£4  +  560£3  -  220£2  +  30£) 

Linear 

Wagner 

0.2£ 

Q  ^ 2  arcsin( V7)  +sin (2  arcsin( V7) ) 

Q  f  sin(2fcarcsin(\/?))  ,  sin(2(fe-l)arcsin(v^))  \  0  „ 

Polynomial 

Rear  loading 

0.52  (0.5£3  —  1.5£2  +  £) 
o.4  (e  -  e) 

0.52  (0.5£  —  0.5£3) 

6625000  (1—0  £15e1/5-2°£ 

17500000  (1-0  ^18e1/4-2°f 

44440000  (1-0  £22.66e7/20-20£ 

90000000  (1-0  ^30e1/2-20^ 

Table  1 :  Modification  functions  used  in  the  design  examples. 


Figure  4:  Optimization  history  and  obtained  shapes  in  the  Euler+Bl  optimization  runs. 
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The  optimization  history,  and  the  initial  and  final  shapes  are  reported  in  figure  4.  The  initial  drag  coefficient 
was  Cd  =  0.009070  at  a  =  1.9216°,  while  the  final  one  resulted  to  be  Cd  =  0.008645  at  a  =  1.5619°.  A 
subsequent  analysis  with  the  ZEN  NS  flow  solver  on  the  same  configurations,  however,  showed  different  trends, 
and  there  was  no  appreciable  difference  between  the  drag  coefficients  of  the  two  configurations.  Therefore,  it 
was  decided  to  use  the  initial  RAE  2822  as  starting  point  for  the  subsequent  optimization  pass  with  the  Navier- 
Stokes  solver. 

5.2  Optimization  using  hybrid  GA  and  NS 

The  result  of  the  previous  run  were  used  to  select  the  design  variables  that  had  the  strongest  effect  on  the  objective 
function. 

A  256  x  56  C-shaped  grid  was  used  for  each  Navier-Stokes  run.  A  three  level  multigrid  acceleration  strategy 
was  used. 


0.01043 


0.01042 

0.01041 

0.01040 

0.01039 

■O 

O 

0.01038 

0.01037 

0.01036 

0.01035 

0.01034 

0  5  10  15  20  25 
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Figure  5:  Optimization  history  with  BFGS  algorithm,  Navier-Stokes  equations  for  drag  evaluation  and  four  de¬ 
sign  variables. 

A  first  set  of  4  design  variables  was  identified  and  a  simple  BFGS  run  was  performed.  The  maximum  num¬ 
ber  of  BFGS  steps  was  fixed  to  10.  The  optimization  history  is  reported  in  figure  5.  The  optimization  process 
terminated  when  the  convergence  criteria  were  meet,  after  three  gradient  evaluations  and  a  total  of  23  objec¬ 
tive  function  evaluations.  As  can  be  observed,  a  small  drag  reduction  was  obtained  (0.7%).  It  is  worth  to  note 
that  in  figure  5,  as  well  as  in  all  the  following  evolution  histories  reported,  all  the  objective  function  evalua¬ 
tions  required  by  the  optimization  algorithm  are  reported,  including  those  used  to  compute  the  gradient  by  finite 
differences. 

The  number  of  design  variables  was  then  extended  to  eight.  Here  three  different  run  were  performed,  The 
first  one  was  a  simple  GA  with  8  population  members  that  ran  for  15  generations.  The  second  one  was  a  hybrid 
GA  that  ran  for  4  generation,  with  a  population  size  of  two  elements.  The  BFGS  was  applied  to  each  element, 
and  the  number  of  BFGS  steps  was  fixed  to  10.  This  produced  a  total  of  62  objective  function  evaluations.  The 
last  hybrid  GA  was  instead  characterized  by  a  population  of  eight  elements  and  the  BFGS  activation  probability 
was  set  to  6%  with  three  descent  steps  allowed.  After  7  generation  the  optimizer  required  82  objective  calcu¬ 
lations.  The  evolution  history  for  these  three  runs  is  reported  in  figure  6.  In  run  3  the  best  solution  had  a  drag 
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coefficient  equal  to  0.010295  while  the  starting  one  was  Cd  =  0.010415.  The  cp  distribution  of  the  original  and 
the  optimized  airfoil  are  reported  in  figure  7,  the  c /  comparison  is  reported  in  figure  8,  and,  finally,  the  Mach 
field  around  the  airfoil  is  reported  in  figure  9. 
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Figure  6:  Evolution  history  of  the  three  runs  with  eight  design  variables  and  Navier-Sokes  solver. 


6  Use  of  the  adjoint  method  for  gradient  computation 

A  planned  improvement  of  the  hybrid  genetic  algorithm  is  the  introduction  of  the  adjoint  method  in  gradient 
computation.  The  major  advantage  of  this  approach  is  that  there  a  higher  computational  efficiency  is  obtained. 
Indeed,  by  solving  the  adjoint  equations,  one  obtains  all  of  the  gradient  components. 

A  multi-block  continuous  Navier-Stokes  adjoint  equation  solver  for  two-dimensional  fields  was  implemented. 
The  numerical  solution  of  the  adjoint  equations  is  obtained  by  using  a  first-order  time-dependent  technique  based 
on  a  finite  volume  discretization.  The  solver  computes  the  fluxes  at  cell  interfaces  using  a  flux-vector  splitting 
technique.  In  a  similar  way,  the  boundary  conditions  are  imposed  on  the  numerical  fluxes  at  the  computational 
field  edges. 

At  the  moment  the  adjoint  based  procedure  is  being  developed  as  a  stand-alone  code,  whose  scheme  is  re¬ 
ported  in  figure  10.  The  gradient  computed  using  the  adjoint  equation  set  is  then  used  by  a  conjugate  gradient 
optimization  routine.  After  the  validation  phase,  the  developed  solver  will  be  used  as  core  of  a  new  hybrid  op¬ 
erator  of  the  genetic  algorithm. 

The  functional  to  be  minimized  is: 

C  =  lo1D  +  lo2(L~^2 

where  D  is  the  drag  L  is  the  lift,  L*  is  the  desired  lift  and  the  ujt  are  weights. 

The  first  case  considereds  corresponds  to  the  the  case  of  fixed  lift  without  constraints  on  the  geometry.  The 
absence  of  constraints  on  shape  will  lead,  as  will  be  evident  in  the  first  example  reported,  to  a  clear  tendency  to 
thickness  reduction  and  to  rear  loading  the  airfoil. 

Figure  11  reports  some  results  related  to  the  design  problem  previously  defined  (RAE2822,  M  =  0.68, 
Re  =  5.7  x  106,  ci  =  0.56)  when  no  control  on  maximum  thickness  is  imposed.  In  particular,  drag,  penalty 
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Figure  7:  Comparison  between  initial  and  final  cp  distributions  in  NS  run  3. 


Figure  8:  Skin  friction  coefficient  comparison  in  NS  run  3. 


on  lift,  and  initial  and  final  airfoil  shapes  are  reported.  The  gradient  is  computed  point-wise,  i.e.  each  grid  point 
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Figure  9:  Mach  field  around  the  optimized  airfoil  of  run  3. 


on  the  airfoil  boundary  is  a  design  variable.  The  gradient  analytic  expression  involves  second  as  well  as  third 
derivatives  of  the  flow  variables.  As  the  flow  field  is  only  second  order  accurate  in  space,  it  can  occur  that  the 
computation  of  such  derivatives  is  not  reliable  in  regions  of  abrupt  flow  changes,  as  for  example  near  the  front 
stagnation  and  close  to  the  rear  separation.  For  this  reason  the  point- wise  gradient  was  smoothed  using  a  Fourier 
filter.  After  such  filtering  the  conjugate  gradient  method  is  able  to  nicely  decrease  the  functional  C. 

Figure  12,  instead,  reports  the  result  obtained  when  an  a-posteriori  control  on  maximum  thickness  is  im¬ 
posed.  This  is  obtained  projecting,  at  each  iteration,  the  airfoil  modification  vector  along  a  direction  in  which 
the  modification  in  the  maximum  thickness  section  is  zero.  It  is  worth  to  note  that  while  in  the  first  run  over  400 
field  computations  were  allowed,  in  this  second  one  only  about  160  evaluations  were  allowed  in  order  to  save 
computational  resources. 


7  Conclusions 


Navier-Stokes  flow  solvers  can  be  profitably  used  for  aerodynamic  shape  design  with  evolutionary  optimizers, 
provided  that  attention  is  paid  to  computational  efficiency.  In  particular,  a  pre-analysis  of  the  design  space  with 
lower  fidelity  tools  is  highly  recommended  to  avoid  waste  of  computational  resources. 

Furthermore,  parallel  computing  and  advanced  approaches  like  gradient  computation  through  adjoint,  pave 
the  way  for  challenging  industrial  applications. 
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Figure  10:  Adjoint  optimizer  scheme. 


Acknowledgments 

Thanks  arc  due  to  prof.  Mark  Drela  for  allowing  the  use  of  his  MSES  code  in  our  research  work,  and  to  dr.  Pietro 
Catalano  and  dr.  Luigi  Paparone  of  the  Italian  Center  for  Aerospace  Research  (CIRA)  for  their  precious  help  on 
the  use  of  ZEN  flow  solver. 

Part  of  the  work  here  presented  was  carried  out  within  the  framework  of  the  AEROSHAPE  project.  This 
project  (Multi-Point  Aerodynamic  Shape  Optimization)  is  a  collaboration  between  Aerospatiale  Matra  Airbus, 
Alenia  Aeronautica  (Coordinator),  Daimler  Crysler  Airbus,  EADS-M,  Dassault  Aviation,  SAAB,  SENER, 
SYNAPS,  CIRA,  DERA,  DLR,  FFA,  INRIA,  HCSA,  NLR,  ONERA,  and  NuTech  Solutions.  The  project  is 
funded  by  the  European  Commission,  DG  Research,  under  the  GROWTH  initiative  (Project  Ref.:  GRD1-1999- 
10752). 

References 

[1]  Shigheru  Obayashi.  Pareto  genetic  algorithm  for  aerodynamic  design  using  the  navier-stokes  equations. 
In  Domenico  Quagliarella,  Jacques  Periaux,  Carlo  Poloni,  and  Gabriel  Winter,  editors.  Genetic  Algo¬ 
rithms  and  Evolution  Strategies  in  Engineering  and  Computer  Science,  pages  245-266,  England,  Novem¬ 
ber  1997.  John  Wiley  &  Sons  Ltd. 

[2]  N.  Marco,  S.  Lanteri,  J.  A.  Desideri,  and  Jacques  Periaux.  A  parallel  genetic  algorithm  for  multi-objective 
optimization  in  computational  fluid  dynamics.  In  Kaisa  Miettinen,  Marko  M.  Makela,  Pekka  Neittaanmaki, 
and  Jacques  Periaux,  editors.  Evolutionary  Algorithms  in  Engineering  and  Computer  Science ,  pages  445- 
455,  England,  1999.  John  Wiley  &  Sons  Ltd. 


1/2(c,-k) 


12-11 


RAE  2822 


Figure  1 1 :  Drag,  lift  penalty  constraint  for  the  functional  C,  and  airfoil  shapes  comparison  for  the  case  without 
constraint  on  maximum  thickness. 
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Figure  12:  Drag,  lift  penalty  constraint  for  the  functional  C,  and  airfoil  shapes  comparison  for  the  case  with 
constrained  maximum  thickness. 
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